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Abstract 

In this second report on our recent numerical simulations of two-flavour QCD, we provide 
further technical details on the simulations and describe the methods we used to extract the 
meson masses and decay constants from the generated ensembles of gauge fields. Among 
the topics covered are the choice of the DD-HMC parameters, the issue of stability, autocor- 
relations and the statistical error analysis. Extensive data tables are included as well as a 
short discussion of the quark-mass dependence in partially quenched QCD, supplementing 
the physics analysis that was presented in the first paper in this series. 



1. Introduction 



Lattice QCD with Wilson quarks [1] has seen important algorithmic developments 
in the last few years [2-8]. As a consequence, a large range of lattice spacings, lattice 
volumes and quark masses can now be explored, using numerical simulations, thus 
providing new physics opportunities and a greater lever arm for the extrapolations 
to the continuum and the chiral limit. Our recent work [9] was the first to fully profit 
from the technical breakthrough and several other projects, simulating QCD with 

* On leave from Centre de Physique Theorique, CNRS Luminy, F-13288 Marseille, France 
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two [10,11] and three [12-14] flavours of light Wilson quarks, or with two flavours 
and a twisted mass term [15], are currently underway, all heavily depending on the 
new generation of algorithms. 

The present paper is the second in a series of two papers devoted to the study of 
two-flavour QCD at small quark masses and lattice spacings. In the first paper [9], 
the focus was on the physics results, while here we give a fairly detailed technical 
account of the simulations that we have performed. 

Perhaps the most important items that we discuss are the stability of the simula- 
tions (sect. 3) and the pattern of autocorrelation times observed in our runs (sect. 4). 
We also describe, in sect. 5, the methods that we used to extract the meson masses 
and decay constants from the generated ensembles of gauge fields (extensive data 
tables are included in appendix C). The paper ends with an addendum to the first 
paper, where we briefly discuss the quark-mass dependence of various quantities in 
partially quenched QCD with 2 + 1 flavours of quarks. 



2. Simulation parameters 

We consider the Wilson formulation of lattice QCD, optionally 0(a)-improved, with 
a doublet of mass-degenerate sea quarks. The notation and normalization conven- 
tions adopted in this paper coincide with those already used in our previous paper 
[9]. In particular, the parameters of the lattice theory are the inverse bare coupling 
f3, the sea-quark hopping parameter K sea _ and the coefficient c sw of the Sheikoleslami- 
Wohlert improvement term [16,17]. 

All simulations reported here were performed using the DD-HMC simulation algo- 
rithm [7]. As suggested by the name, the algorithm combines domain-decomposition 
ideas with the HMC algorithm [18]. More precisely, by dividing the lattice into non- 
overlapping rectangular blocks, a natural separation of the high-frequency from the 
low-frequency modes of the fields is achieved. Following Sexton and Weingarten [19], 
the different modes are then evolved using different molecular-dynamics step sizes, 
which results in a significant acceleration of the simulation. 

On a given lattice and at fixed coupling, the simulations progressed from the larger 
to the smaller quark masses, normally skipping 1500 molecular-dynamics trajectories 
for thermalization. The number iV tr j of trajectories generated after thermalization, 
the separation N sep (in numbers of trajectories) between successive saved field con- 
figurations and the number -/V c f g of saved fields are given in table 1. Different runs at 
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Table 1. Lattice parameters and simulation statistics 



Run 


Lattice 


P 




^sea 




iVsep 


iVcfg 


A 

Ma 


32 X 24 


0.0 


(J 


0.15750 


conn 

6300 


1 nn 

100 


a a 


A 








0.15750 


5070 


on 

30 


1 c n 

169 


A 

M 








0.15800 


1 nonn 
10800 


1 nn 
100 


1 nn 
109 


A 

A 3 a 








0.15825 


6100 


100 


62 


A 

A 3 b 








0.15825 


3800 


1 nn 
100 


no 

38 


A 

A A 








0.15835 


4950 


£n 
50 


1 nn 
100 


r> 
Dl 


64 X o2 


O.O 


u 


U.1541U 


5050 


£n 
50 


1 nn 
100 










n 1 c/i /in 
U.1544U 


rnnn 

5200 


£n 
50 


1 ni 
101 


B 3 








0.15455 


5150 


50 


104 










0.15462 


5050 


50 


102 


c\ 


64 x 24 3 


5.6 





0.15800 


3450 


30 


116 


Di 


48 x 24 3 


5.3 


1.90952 


0.13550 


5150 


50 


104 


D 2 








0.13590 


5130 


30 


171 


D 3 








0.13610 


5040 


30 


168 


D 4 








0.13620 


5010 


30 


168 


D 5 








0.13625 


5040 


30 


169 


Ei 


64 x 32 3 


5.3 


1.90952 


0.13550 


5344 


32 


168 


E2 








0.13590 


5024 


32 


158 


E 3 








0.13605 


5024 


32 


158 



the same lattice parameters (such as A 3a and A 3o ) are distinguished by a lower-case 
latin index. In our previous paper [9], only the runs A\ a , A 2 , A 3a , A 3 b, B1-B4 and 
Di~D 5 were included in the physics analysis. The other runs listed in table 1 merely 
serve, in sections 3 and 4, to clarify some technical issues. 

The DD-HMC simulation algorithm was implemented following the lines of ref. [7] . 
In particular, for the solution of the Dirac equation on the full lattice, the Schwarz- 
preconditioned GCR solver described in ref. [6] was used. The so-called replay trick, 
however, was switched off in the more recent simulations A 3 b~E 3 , because trajectory 
replays would have been rare and hardly worth the extra effort (see subsect. 3.3). 

No attempt was made to tune the DD-HMC algorithm and most of its parameters 
were actually set to some fixed values, the same as the ones already chosen in ref. [7]. 
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Table 2. DD-HMC parameters, acceptance rate and average solver iteration numbers 



Run 


Block size 


N 2 


p 

1 acc 


(^gcr) 


(A^cg) 


Ai 


8 x fi 2 x 1 2 




81 * 


23 


73 


■"■lb 


8 x 6 x 12 2 


5 


0.82 


22 


89 




8 x 6 2 x 12 


6 


0.79* 


39 


74 






10 


0.89* 


54 


75 






10 


0.86 


54 


75 


Aa 




1 fi 


91 


73 


75 

1 'J 




8 3 x Ifi 


s 

O 


84 


32 


85 


J->2 




1 


89 


52 


87 


B 3 




12 


0.87 


74 


87 


B A 




14 


0.92 


90 


88 


Ci 


8 x 6 x 12 2 


7 


0.81 


41 


92 


Di 


6 2 x 12 2 


7 


0.81 


25 


120 


D 2 




8 


0.80 


41 


123 


D 3 




12 


0.87 


58 


124 


D 4 




14 


0.87 


73 


125 


D 5 




18 


0.89 


87 


125 


E 1 


8 4 


9 


0.80 


25 


121 


E2 




11 


0.84 


41 


124 


E 3 




13 


0.83 


53 


125 



* Transition probability includes trajectory replays 



Among these were the trajectory length r = 0.5, the integration step numbers 
Nq = A and iVi = 5 associated to the gauge and block fermion forces as well as the 
admitted tolerances (n, r%, f\, = (10 -8 , 10 -7 , 10 -11 , 10 -10 ) for the numerical 
solution of the Dirac equation on the blocks and the full lattice f. The parameters of 
the Schwarz-preconditioned GCR solver were fixed to the values quoted in ref. [6], 

f The trajectory length r and thus the integration step sizes t/N^, etc., refer to a particular nor- 
malization of the kinetic term in the molecular-dynamics Hamiltonian. Here the normalizations are 
the same as in ref. [7], i.e. the term is assumed to be equal to i(n,II) = tr{II(x, ii)^Il(x, /i)}, 

where Yl(x,i\j,) denotes the canonical momentum of the link variable U(x,fi). 
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except for the number n^v of Krylov vectors generated before the GCR recursion is 
restarted, which was set to 32 in run D§ and to 24 in all other runs. 

What remains to be specified are then the size of the blocks on which the algorithm 
operates and the integration step number N2 associated to the block interaction term 
in the molecular-dynamics Hamiltonian (see table 2). In practice the latter must be 
increased as one moves to lighter quark masses in order to preserve a high acceptance 
rate P acc . The average number Nqqr of GCR solver iterations needed along the 
trajectories also depends on N2 (it decreases when N2 goes up), while the average 
number Nqg of conjugate-gradient iterations required for the computation of the 
block terms in the molecular-dynamics equations is largely determined by the block 
size. 

With the chosen parameters, the reversibility of the molecular-dynamics trajecto- 
ries is guaranteed to high precision. In the tests that we have performed, the average 
absolute deviation of the components of the link variables after a return trajectory 
was at most 3 x 10~ 9 , while in the case of the Hamiltonian the observed differences 
were less than 4 x 10 -6 . Deviations larger than 10 times the average occurred in 
less than 1% of the cases and never went beyond 100 times the average. 



3. Spectral gap and stability issues 

The Wilson-Dirac operator preserves chiral symmetry only up to lattice effects and 
is therefore not rigorously protected from having eigenvalues much smaller than the 
quark mass. Exceptionally small eigenvalues do not invalidate the theory but may 
lead to instabilities in numerical simulations, depending, to some extent, on which 
simulation algorithm is used. 

3. 1 Spectral gap of the Dirac operator 

In a previous dedicated study [20] , we computed the distribution of the spectral gap 
of the hermitian lattice Dirac operator on the lattices A\ — A4, B\, B2, C\ and D\. 
The distributions turned out to be well separated from the origin, thus showing, a 
posteriori, that the simulations were safe of exceptionally small eigenvalues and the 
associated instabilities. Moreover, based on the observed scaling properties of the 
distributions on the A, B and C lattices, we argued that this will always be so in 
the large-volume regime of the unimproved Wilson theory. 

The gap distributions have now also been computed on the lattices B3, B4, D2 — 
D 5 , E 2 and E 3 . In the following, however, we focus on the improved theory, because 
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Fig. 1. Normalized histograms of the (unrenormalized) spectral gap fi of the her- 
mitian lattice Dirac operator, as obtained in the runs D2 — D5 . The bin size is 2 MeV 
and the dotted vertical lines indicate the position of the median p, of the distributions. 
The data were converted to physical units using a — 0.0784 fm [9]. 



the results obtained on the B lattices are fully in line with the behaviour expected 
from our previous paper [20]. 

At first sight, the gap distributions in the improved theory look similar to those 
in the unimproved theory (see fig. 1). In particular, they are well separated from the 
origin, on all lattices that we have simulated, and the median of the distributions 
again turns out to be a practically linear function of the sea-quark mass (fig. 2). 

However, the dependence of the width a of the distributions on the quark mass 
and the lattice size is different (see table 3) f . In the case of the D-series of lattices, 
for example, the width decreases by as much as a factor of 1.5 from the largest to 
the smallest quark mass, while no obvious mass-dependence was seen on the A and 

I Following ref. [20], we define the width of the distributions through a = h(v — u), where [u,v] 
is the smallest range in \l, which contains more than 68.3% of the data. 
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Fig. 2. Median p, of the gap distributions obtained in runs D\ — D 5 (data points), 
plotted as a function of the bare sea-quark mass m sea (see subsect. 5.2 for the precise 
definition of the latter) . The line is a linear fit of the data without constant term. 

B lattices. Moreover, a does not appear to scale proportionally to the inverse square 
root of the (four-dimensional) volume V of the lattice (see fig. 3). The widths on 
the lattices Di and E2, for example, turned out to be nearly the same, contrary to 
what was expected on the basis of the experience made in the unimproved theory. 

Another perhaps not unrelated observation is that the median of the distribution 
on the D and E lattices is always smaller than the threshold of the spectral density 
in infinite volume, which we expect to be at Z^m sea _ [20], Zj± being the axial-current 
renormalization constant (Z& = 0.75(1) on these lattices [21]). The spectral density 
in finite volume thus has a tail that extends a few MeV below the threshold. On the 
other hand, the values quoted in table 3 of the average splitting (A) of the lowest 
four eigenvalues suggest that the tail scales to zero in the infinite- volume limit, as it 
has to be if the density in infinite volume does not extend all the way to zero [20]. 

At present, however, there is still no theoretical understanding of the dependence 
of the gap distribution on the quark mass and the lattice size. In particular, the fact 
that the improved and the unimproved theory behave differently in this respect re- 
mains unexplained. Partially quenched (Wilson) chiral perturbation theory may be 
a framework in which these questions can be addressed [22] and further insight may 
perhaps also be gained by studying the localization properties of the eigenfunctions 
and the convergence of the spectral density to the infinite-volume limit. It would 
be interesting to know, for example, whether the spectral gap coincides with the 
mobility edge [23] and whether the tail of the spectral density below Zpjn 9 ^ does 
in fact disappear in the infinite- volume limit. 
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Table 3. Median and width of the gap distributions in the improved theory 



Run 



a 



IX - Z A m. 



'sea 



(A) 



D 2 

D 5 
E2 
E 3 



57.3(6) 
32.0(3) 
21.4(3) 
15.9(3) 
12.9(4) 
30.3(3) 
21.3(3) 



3.3(4) 

2.79(24) 

2.84(23) 

2.33(18) 

1.99(15) 

2.58(19) 

2.31(19) 



6.5(10) 

4.8(6) 

2.3(4) 

2.1(3) 

1.4(4) 

6.6(6) 

5.8(5) 



2.48(12) 

2.39(7) 

2.29(7) 

2.23(6) 

2.28(5) 

1.69(8) 

1.52(7) 



All entries are given in MeV 



3.2 Accessible range of pion masses on the D and E lattices 

When the sea-quark mass decreases, the gap distribution becomes sharper and moves 
closer to the origin. Eventually the probability for exceptionally small eigenvalues 
is not completely negligible anymore and one may run into algorithmic instabilities. 
We have not reached this point yet and consequently cannot say in which way the 
DD-HMC simulations will be affected. However, in order to be on the safe side, 
one may prefer to stay in the range of parameters where the gap distribution is well 
separated from the origin, i.e. where, say, the inequality fx > 3a holds [20]. 

On a given lattice, this bound sets a lower limit on the accessible sea-quark masses 
and thus on the masses M w of the pions (the lightest pseudo-scalar mesons made of 
the sea quarks). Furthermore, if large finite- volume effects are to be avoided, the 
bound M n L > 3 (where L denotes the spatial lattice size) should better be respected 
as well. 

In the case of the D and E lattices, the range of pion masses where both conditions 
are fulfilled can be determined explicitly, using our simulation results. An extra- 
polation in the sea-quark mass is however still required, but it seems reasonable to 
extrapolate jx and M% linearly [9] and to assume that a drops to values below 2 MeV 
at small quark masses. For the accessible range of pion masses, we then obtain 



where the limit is set by the constraint M^L > 3 on the D lattices. This is not so on 
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Fig. 3. Width a of the gap distributions, scaled by the factor y/V /a, as obtained 
in the unimproved (left) and the improved theory (right). The statistical errors were 
determined using the bootstrap method. 

the E lattices, but values of M^L as low as 3.4 can still be safely reached, i.e. also 
in this case, the stability bound is not too restrictive. 



3.3 Molecular- dynamics instabilities 

Similar to the standard HMC algorithm, the DD-HMC algorithm obtains the next 
field configuration by integrating the associated molecular-dynamics equations. The 
numerical integration of these equations is well known to be potentially unstable. 
If an instability occurs, the energy deficit AH at the end of the integration can be 
large and the new field configuration is then normally rejected. The efficiency of the 
simulation may thus be affected, but we wish to emphasize that large energy deficits 
do not invalidate the algorithm unless the reversibility of the molecular-dynamics 
integration is compromised. 

Earlier studies of the phenomenon suggest that the instabilities are caused by ex- 
ceptionally small eigenvalues of the lattice Dirac operator [24-26] . Even if the gap 
distribution is safely separated from zero, it is possible that the Dirac operator de- 
velops such eigenvalues somewhere along the molecular-dynamics trajectories. The 
probability for this depends on how accurately the molecular-dynamics equations 
are solved, i.e. on the integration step sizes and the solver residues. 

In our simulations, the probability for \AH\ to be larger than 2 was always fairly 
small and often equal to zero (runs Bi — B4, for example). The worst cases in the 
unimproved and the improved theory were the runs A4 and D 5 respectively, where 
the threshold of 2 was passed by 1.4% and 0.7% of the trajectories. Energy deficits 
\AH\ larger than 10 3 were never seen, but values above 100 did occur, although very 
rarely so. 
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4. Autocorrelation times 



The dynamical properties of the simulation algorithms used in lattice QCD are still 
largely unknown. It is not clear, for example, whether there are several relevant time 
scales and how they depend on the lattice parameters and the chosen algorithm. We 
shall not attempt to answer these difficult questions here, however, and merely give 
an account of our empirical studies of the autocorrelations in the runs listed in 
table 1. 

4-1 Determination of autocorrelation times 

Following the standard conventions, we define the integrated autocorrelation time 
Tint of an infinite series 01, 02, 03, . . . of measured values of an observable A through 

where T(t) denotes the autocorrelation function of the series. In practice only a finite 
number N of measurements can be made and the estimation of the autocorrelation 
time from the available data then requires some ad hoc choices to be made. 
For the autocorrelation function we use the approximation 

j N-t 

r ^ ~ W~t E ( ai ~ - «+)> < t < iV, (4.2) 

i=l 

in which a_ and a+ are, respectively, the averages of the first N — t and the last 
N — t elements of the series ai, . . . , a at. The sum in eq. (4.1) is then truncated at 
some value W <C N of the time lag t, referred to as the summation window, which 
should ideally be such that the remainder of the sum can be safely neglected. 

If the autocorrelation function is well behaved, as in the case shown in the upper 
plot of fig. 4, the choice of the summation window is not critical and any reasonable 
prescription will do. The rule adopted here is to stop the summation in eq. (4.1) 
at the first value of t where the normalized autocorrelation function is equal to zero 
within two times its statistical error, the latter being estimated using the Madras- 
Sokal approximation (see appendix E of ref. [7]). 

In practice the calculated autocorrelation functions may have long tails and they 
may also vary significantly with the selected range of the data series. An example 
illustrating this behaviour is shown in the lower plot in fig. 4. In all these cases, we 
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Fig. 4. Normalized autocorrelation functions T(t)/T(0), plotted versus the time 
lag t given in numbers of trajectories, of the plaquette P (upper plot) and the solver 
iteration number Nqcr (lower plot). The data shown were calculated using the last 
4000 trajectories of run B2 (full points) or only the first 2000 of these (open points). 

divide the data series into large bins, calculate the bin averages and estimate the sta- 
tistical variance a 2 of the total average assuming these are statistically independent. 
The integrated autocorrelation time is then given by 



Tint 



2aV 



(4.3) 



where do denotes the naive statistical error. Evidently, the results obtained in this 
way are rough estimates that could easily be wrong by factor 2 or so. 

4-2 Reference autocorrelation times 

The integrated autocorrelation times of the Wilson plaquette P and the GCR solver 
iteration number Nqqr are listed in table 4. These two quantities are unphysical, 
but they are readily accessible and are useful reference cases that probe the dynamics 
of the simulation at both short and long distances. 

In order to facilitate the comparison of the figures quoted in the table, the auto- 
correlation times were determined using data series of a fixed length equal to 4000 
trajectories. The autocorrelation times are given in numbers of trajectories and er- 
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Table 4. Autocorrelation times of the plaquette P and the solver iteration number Aqcr 



Run 


Ti„t[P] 


T int [^GCr] 


Run 


Ti„t[P] 


Tint [-NGCR,] 


Ma 


25(5) 


43* 


Ci 


17(3) 


35(7) 


A lb 


29(6) 


38* 


D 1 


11(1) 


10(2) 


A 2 


23(4) 


46* 


D 2 


17(3) 


21(4) 


A 3a 


14(2) 


53(10) 


D 3 


16(2) 


19(3) 


A 3b 


28* 


53* 


D 4 


16(2) 


15(2) 


A A 


19(4) 


45* 


D 5 


32(6) 


24(5) 


Bi 


14(2) 


50* 


E 1 


33* 


14(3) 


B 2 


12(2) 


39* 


E 2 


19(3) 


11(2) 


B 3 


9(1) 


45* 


E 3 


27(5) 


25(5) 


B 4 


14(2) 


51* 








* Estimate based on 


data binning 









ror estimates are quoted only in those cases where the autocorrelation function was 
well behaved. In these regular situations, the binning method always gave consistent 
results. 

In all simulations of the improved theory, except for run E\ perhaps, the autocor- 
relation times were safely determined and turned out to be reasonably small. This 
was not so in the simulations of the unimproved theory, where the autocorrelation 
function of the GCR iteration number typically had a tail similar to the one shown 
in the lower plot in fig. 4. 0(a) improvement thus appears to have the side-effect of 
reducing the autocorrelation times. 

The regularity of run C\ then remains unexplained, however, and the differences 
in the autocorrelation times could actually also very well be related to the fact that 
the physical volumes of the C\ — E 3 lattices are larger, by a factor of two or more, 
than the volumes of the other lattices. Presumably the size of the blocks, on which 
the DD-HMC algorithm operates, matters as well, although the comparison of the 
runs Ai a and An, does not suggest this to be so. 

4-3 Autocorrelations of physical quantities 

The meson masses and all other physical quantities were calculated after finishing 
the simulations, using the generated ensembles of saved gauge-field configurations 
(see table 1). A fairly large number of trajectories was skipped between successive 
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saved configurations so that the statistical correlations in these sets of fields can be 
expected to be small. 

In order to find out whether the residual correlations are relevant for the de- 
termination of the statistical errors, the basic two-point correlation functions were 
averaged over small bins of successive configurations. The physical quantities were 
then extracted from the binned data and their statistical errors were estimated using 
the jackknife method (appendix A). If there were significant statistical correlations 
in the data, the errors would increase with the bin size, but this was not the case 
and we therefore concluded that it was safe to proceed without data binning. 



5. Computation of meson masses and decay constants 

The masses and matrix elements tabulated in appendix C were calculated using a 
combination of methods, most of which being entirely standard by now. We consider 
two valence quarks, labelled r and s, and study the vector and pseudo-scalar mesons 
in the fs-channel. The masses of the valence quarks may be set to the sea-quark 
mass, but we are also interested in the partially quenched situation where one of the 
quark masses is different from the sea-quark mass. 

5. 1 Two-point correlation functions 

The pseudo-scalar density, the axial current and the vector current in the fs-channel 
are given by 

P rs = n 5 s, A™ = f 7^755, V™ = fj„s. (5.1) 

All masses and decay constants we are interested in were extracted from the two- 
point functions 

/pp(x )=« 3 £ (P rs (x)P sr (0)), (5.2) 

x lt x 2 ,x 3 

/ A p(x ) = a 3 (A r o s (x)P sr (0)), (5.3) 

Xi,X 2 ,X 3 

3 

hv(xo)=a 3 J2 J2(W r k s (x)W^(0)), (5.4) 

x±,x 2 ,x 3 k=l 
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where W™ is a linear combination of the vector current VT S and a Jacobi smeared 
form of it [27] , slightly tuned so as to suppress the high-energy intermediate states 
in the two-point function. 

The correlation functions were evaluated in the standard manner by first express- 
ing them as an expectation value of a product of two quark propagators. These were 
calculated by solving the lattice Dirac equation, using the Schwarz-preconditioned 
GCR solver [6] and requiring the normalized residue of the solution to be less than 
10 -10 . In order to reduce the statistical fluctuations, the results were averaged over 
time-reflections and 5 distant source points in the case of the A and B runs and over 
3 source points in the case of the D runs. 

5.2 Masses and matrix elements 

On a lattice of infinite time- like extent, and at large times xq, the correlation function 
fpp(xo) is saturated by the one-particle pseudo-scalar meson state in the fs-channel. 
If we denote the mass of the meson by Mps and the associated vacuum-to-meson 
matrix element by Gps, the asymptotic form of the correlation function is 

/pp(x ) = -^e- M — « + ..., (5.5) 

where the ellipsis stands for a series of more rapidly decaying terms. The mass My 
of the fs vector meson may be defined similarly through the asymptotic behaviour 
of the vector correlation function fw(xo), but the definition requires further expla- 
nation if the meson is unstable in infinite volume (see subsect. 5.6). 
Next we note that the ratio 

m eS (x ) = {\ (d + do) /ap(zo) + CAa<9o<9 /pp(xo)} / fpp(x ) (5.6) 

converges to a constant m rs at large times xo, for any fixed value of the parameter 
ca, because both /ap(xo) and fpp(xo) are proportional to e _Mpsa::o in this limit. 
Moreover, in the continuum limit, m e ff(xo) is expected to converge to the sum of 
the bare current-quark masses of the r and the s quark, at all times x , with a rate 
proportional to a in the unimproved theory (where we set ca to zero) or a 2 if the 
improvement coefficients c sw and ca are properly tuned [17,28,29] f. 

All our numerical data for m e s(xo) in fact turned out to be statistically consis- 
tent with a constant value, over a large range of x , and the quark mass sum m rs 

| The effects of the 1 + O(am) renormalization factors (C.2) are expected to be small in practice 
and are neglected here for simplicity. 
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was therefore always unambiguously and accurately determined. In particular, the 
current-quark of the sea quarks is obtained by setting the hop- 

ping parameters of the valence quarks to K sea . Whether in general m rs coincides 
with \{m rr + m ss ), as one expects to be the case if the lattice effects are small, is 
a question to which we shall return in sect. 6. 

The bare pseudo-scalar decay constant Fp$ in the fs-channel is normally extracted 
from the asymptotic behaviour of the two-point functions /ap(^o) an d /pp(a>o)- in 
this paper, however, we first computed m rs , Mps and Gps and then used the formula 

F PS = |^Gps (5.7) 

for the decay constant. Starting from eq. (5.6), it is straightforward to show that 
equivalent results are obtained in this way, up to small corrections of 0(a 2 ). Note 
that Fps is automatically 0(a)-improved if m rs is. 

5. 3 Spectral decomposition in finite volume 

On a finite lattice with time-like extent T, the calculation of the pseudo-scalar and 
vector meson masses requires some care and must address the issue of higher-states 
contributions. This is, incidentally, not so in the case of the quark mass sum m rs , 
which is expected to be independent of the lattice size up to lattice-spacing effects. 

For < xo < T, the correlation function /pp(xo) (and similarly /vv(^o)) can be 
expanded in a rapidly convergent series of the form 

oo oo 

fpp(x ) = -^^Cijhixo; Ei, Ej), (5.8) 

i=0 j=i 

h(t; E, E') = exp{-Et - E'(T - t)} + exp{-E't - E(T - t)}, (5.9) 

where = Eq < E\ < E 2 < ■ ■ ■ are the intermediate-state energies and Cjj > the 
associated spectral weights J. In the channel considered here, the lowest intermediate 
state is the fs pseudo-scalar meson state at zero spatial momentum. Then come the 
multi-meson scattering states and more and more complicated states as one moves 
up the energy scale. 

| Equation (5.8) assumes the existence of a positive hermitian transfer matrix which may not be 
guaranteed in the improved theory. It seems likely to us, however, that a transfer matrix can still 
be defined, as is the case in 0(a 2 )-improved gauge theories [30], although complex energy values 
and negative weights may occur at energies on the order of the cutoff scale 1/a. 
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At large xq and T, the dominant term in the series (5.8) is thus the one where Ei = 
and Ej = Mps- Moreover, using the product inequality (B.3), the contributions of 
all higher-energy states can be shown to be exponentially suppressed with respect to 
this term. In practice their effects are seen in the simulation data only when either 
xo or T — xq is not too large. The leading terms in this range are then 

fpp(x ) =c h{x o ;0,M o ) + c 1 h{x o ;0,M 1 ) + M = M PS , (5.10) 

where Mi denotes the energy of the next-to-lowest state in the fs-channel (if the 
spatial volume of the lattice is large enough, this will be a three-meson state with 
all particles at rest). 

Note that each term in the spectral series (5.8) decreases exponentially in the 
range < xq <C \T, with an exponent equal to Ej — Ei that can be as small as the 
pseudo-scalar meson mass, for example, even if both Ei and Ej are not small. The 
presence of such contributions complicates the analysis of the correlation functions 
considerably unless the time-like extent T of the lattice is sufficiently large to strongly 
suppress them. This condition was barely satisfied in the case of the run A4, which 
is why we decided to discard it from the physics analysis (as already mentioned in 
the first paper in this series). 

5.4 Effective masses and matrix elements 

Slightly departing from what is usually done, we define the effective pseudo-scalar 
meson mass M e s(xo) in the fs-channel to be the value of M > where 

hjxp - q;0,M) _ fpp(x - a) 
h{x ;0,M) ~ /pp(x ) 

Using the results obtained in appendix B, it is not difficult to prove that this equation 
has one and only one solution. Moreover, with this definition of the effective mass 
it is guaranteed that M e ff(xo) = Mps at large xq, up to exponentially small terms. 
We then also introduce the effective matrix element 

which converges to Gps in the large-time limit. 

The asymptotic behaviour of the effective mass at large xq and T can be worked 
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(5.11) 



(5.12) 



out explicitly, starting from the spectral representation (5.10). Setting 



^ KSlw ^o> = {^i°%o;0,m)} w (5.13) 

and going through a few lines of algebra, it is straightforward to derive the expansion 

<5 ' 14) 

where the ellipsis stands for terms that are exponentially small with respect to the 
next-to-leading term. A similar formula, 

G eff (so) = G PS |l + ie(x ) + § (1 - *(*(>)) ^jl^° -a) + ' ' ] ' (5 ' 15) 
is obtained in the case of the effective matrix element. 
5. 5 Fit procedures 

From the point of view of the statistical error analysis, the correlation functions /pp, 
/ap and /w are the primary quantities, while the effective quark mass sums, meson 
masses and matrix elements are functions of these. The statistical errors of all these 
quantities tend to be strongly correlated. We took the correlations fully into account, 
from the primary quantities to the final results, by propagating the errors using the 
jackknife method (appendix A). In particular, fitted and interpolated values were 
always considered to be functions of the input data, which allows their errors to be 
calculated in the standard manner. 

The quark mass sum, the pseudo-scalar meson masses and matrix elements, and 
the masses of the vector mesons were all obtained by fitting the corresponding ef- 
fective quantity P e s(xo) in a range to < xo < ii of time with the chosen fit function 
3>(xo). We performed correlated least-squares fits, where the values of the fit pa- 
rameters were determined by minimizing 

ti 

x 2 = E [^effOco) - $(x )] (C" 1 )^ [PefrM - Hvo)] , (5.16) 

xo,Vo=to 

the matrix C being the statistical error covariance of P e f[(to), . . . , P e fr(*i)- The quark 
mass sum m rs , for example, was computed by fitting m e s {xq) to a constant as shown 
in fig. 5a. 
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Fig. 5. Sample plots illustrating the dependence on Xo/a of the effective quark mass 
sum (figure a), the pseudo-scalar mass and matrix element (figures b and c) and the 
vector meson mass (figure d), all given in lattice units. The data points shown are 
from run D4 and the valence quark masses were both set to the sea-quark mass in 
this example. The solid lines are the fits discussed in the text. 



In the case of the pseudo-scalar meson masses, we fitted the data with the asymp- 
totic expression (5.14). We first calculated the mass M n of the pions, i.e. the mesons 
made of the sea-quarks, by substituting Mi = 3M n for the energy of the next-higher 
state (thus assuming the latter is a three-pion state with small interaction energy) 
and adjusting M n and c\/cq so as to minimize x 2 - While the fit curves obtained in 
this way represent the data very well, it should be noted that the fitted value of M T 
is largely determined by the data at large times xq, where a fit to a constant would 
give nearly the same results (see fig. 5b). 

Once M n was determined, the mesons made of a sea quark and a valence quark 
with a mass different from the sea quark were considered. Here we set Mi = Mpg + 
2M n and otherwise proceeded as in the degenerate case. Next the matrix elements 
Gps were computed by fitting the data with the asymptotic expression (5.15), using 
the same values of Mi as in the fits of the effective meson masses (fig. 5c) . We did 
not set ci/cq to the previously computed values, but it turned out that the two fits 
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gave consistent results for this parameter. 
5.6 Energy spectrum in the vector channel 

At small sea-quark masses, the vector mesons become resonances that decay into two 
(or more) pseudo-scalar mesons. As was shown long ago [31], resonances give rise 
to a characteristic volume-dependent pattern of the energy spectrum which allows 
their masses and decay widths to be determined, in principle, from simulation data. 

As before, we considered the channels where one or both of the r and s quarks 
is a sea quark. Starting from the correlation functions fvv(xo), the lowest energy 
My in this channel was calculated by fitting the effective mass with the asymptotic 
formula (5.14) (with Mps replaced by My). For the lowest excited-state energy we 
substituted 

Mi = (M 2 s + k 2 y/ 2 + (M 2 + k 2 ) 1 / 2 , k = 2ir/L, (5.17) 

in this case, L being the spatial size of the lattice. Excellent fits were obtained with 
this ansatz and My was determined quite accurately on all lattices. 

We refer to the energy values My as the vector meson masses in this paper, even 
in those cases where the meson is likely to become a resonance in the infinite volume 
limit (we estimate this to be so at the lightest quark masses in each series of lattices 
and perhaps at some of the second-to-lightest as well) . This use of language is only 
slightly incorrect, however, because in all our simulations My turned out to be at 
most 20% larger than Mps + M n and significantly smaller than Mi , in which case 
the true resonance energy is expected to be close to My [31]. 

We finally note that the statistical errors in the vector channel tend to be larger 
than those in the pseudo-scalar channel. The effect could be related to the resonance 
character of the vector mesons and it is conceivable that a coupled channel analysis, 
such as the one recently presented by Aoki et al. [32] , will not only allow the vector 
meson decays to be studied but may also help to reduce the statistical errors. 



6. Quark-mass dependence in partially quenched QCD 

The most important physical results of our simulations were already presented in 
our first paper in this series [9]. We now discuss the dependence of the quantities 
tabulated in appendix C on the quark masses in some further detail, focusing on the 
empirical facts rather than on their possible theoretical interpretation. 
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Fig. 6. Results for the ratio R n and the difference R^ — Rpg obtained on the D- 
series of lattices. The solid lines represent the global linear fit (6.2). Note that the 
points in the lower plot do not have to line up within errors, since R n — Rvs is a 
function of two independent variables rather than of m se a — ?Ti. va i alone. 



As before we set m sea = ^m rr if the r quark is a sea quark and we now also set 
m vai = ^m S s if the s quark is a valence quark. The figures in the tables are all for 
the mixed case, where one quark is a sea quark and the other a valence quark. We 
are thus considering partially quenched QCD with 2 + 1 flavours of quarks. 

6.1 Quark and pseudo-scalar meson masses 

We first remark that the quark mass sum m rs turns out to be equal to m sea + m va i 
within statistical errors, on all lattices and for all quark-mass combinations. The 
ratio rn rs /(m se!i + m va i) is obtained with better statistical precision than the quark 
masses, but the largest deviation seen in this case is only 0.6%. The additivity of 
the current quark masses (which is an exact property of the theory in the continuum 
limit) is thus accurately guaranteed on the lattices that we have simulated. 
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Fig. 7. Dependence of the bare pion decay constant F n = -Fps^ 1=m and of the 
difference — Fps on the quark masses, as determined on the D-series of lattices. 
The solid lines represent the global linear fit (6.3). 



Next we consider the ratios 

^PS — ; 5 U 7T 



^sea H~ ^val 



R 



Ml 



PS 



2m s 



(6.1) 



which are independent of the quark masses to lowest order of chiral perturbation the- 
ory. However, this is not so at next-to-leading order and the numerically calculated 
ratios are in fact weakly mass-dependent (see fig. 6). An empirical fit 



^Rps = a + ai(m sea + m va i) + a 2 m s 



(6.2) 



represents the data quite well in the given range of masses except perhaps for the 
points where m va i m sea . In the case of the D-series of lattices, for example, the 
data for Rps deviate from the fit by no more than 2% and most points are within a 
margin of 1%. 
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6.2 Pseudo-scalar decay constant and vector meson mass 

As can be seen from the tables in appendix C, the calculated values of Fps/My are 
nearly independent of the quark masses. This comes a bit as a surprise, and could 
merely be an accidental agreement in a limited range of masses, since there does 
not appear to be any obvious physical connection between the pseudo-scalar decay 
constant and the vector meson mass. 

The mass dependence of these two quantities is thus practically the same and it 
suffices to consider one of them. Focusing on the decay constant, a simple linear 
expression, 

Fps = b + b 1 (m sea + m va i) + b 2 m se£l , (6.3) 

turns out to fit the available data for Fp$ very well. On the D-series of lattices, for 
example, the fit matches the data within statistical errors and the maximal relative 
deviation in the given range of masses is only 1.6% (see fig. 7). 

It is tempting to use these fits to extrapolate the decay constant to the chiral limit, 
but as already emphasized in our previous paper [9], such extrapolations are difficult 
to justify and asymptotically inconsistent with one-loop chiral perturbation theory. 
On the other hand, the observed linearity of the pseudo-scalar decay constant in 
the range of masses covered by the simulations is striking and calls for a theoretical 
explanation. 



7. Concluding remarks 

Numerical lattice QCD is currently in an interesting transition phase. The valence 
approximation is now practically overcome, but important physical effects of the 
light sea quarks, such as the decay of the rho meson or the anomaly-driven mass 
splitting between the eta and the pions, still have not or only barely been studied 
directly. Simulations at smaller quark masses and on larger lattices than reported 
here will probably be required for this. Our experience however suggests that the 
prospects for such simulations, using 0(a)-improved Wilson quarks, are now quite 
good. 

So far the DD-HMC algorithm performed well and we did not run into any instabil- 
ities or other technical difficulties. As one moves to smaller quark masses and smaller 
lattice spacings, there may be some room for further algorithmic improvements, but 
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the development of variance-reduction methods is likely to be more rewarding at 
this point, particularly so if disconnected quark-line diagrams and multi-particle 
amplitudes are to be computed. 

The numerical simulations were performed on PC clusters at CERN, the Centra 
Enrico Fermi, the Institut fur Theoretische Physik der Universitat Bern (with a con- 
tribution from the Schweizerischer Nationalfonds) and on a CRAY XT3 at the Swiss 
National Supercomputing Centre (CSCS). We are grateful to all these institutions 
for the continuous support given to this project. 



Appendix A. Statistical error analysis 

In the physics analysis of the runs Ai — A3, B\ — and D\ — D§, we kept track 
of the statistical errors using the jackknife method. In particular, any correlations 
among the errors of different observables were always properly taken into account. 
Here we summarize our conventions and briefly explain the basic procedures that 
we used. 

A.l Jackknife samples 

Let A r , r = 1, . . . , R, be a set of primary stochastic observables and o rj i, • • • , ar,N & 
sequence of N measured values of these. In lattice QCD the most common primary 
observables are the Wilson loops and sums of products of quark propagators. The 
jackknife method assumes that the measured values are unbiased and statistically 
independent. We shall thus take it for granted that the residual autocorrelations are 
negligible in the cases of interest (see sect. 4). 

The averages a r of the observables A r and the associated statistical error covari- 
ance C rs are given by 




(A.l) 



rs 



N(N-l) 



1 




(A.2) 



i=l 



If we introduce the jackknife samples 



a: 



J 

'r,i 



a r + ctv (av — a r,i) ) 



c N = (N(N-l)) 



-1/2 



(A.3) 
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an equivalent expression for the error matrix is 



N 



C rs = ^ ( Q r,i ~ Q r) ( a li ~ a s) ■ (A.4) 
i=l 

Note that our definition of the jackknife samples slightly departs from the standard 
conventions, where cn = 1/(JV — 1). The modification is numerically insignificant 
in practice, but leads to some simplifications when data from different simulations 
are to be combined (see subsect. A. 3). 

A . 2 Error propagation 

Apart from estimating the primary observables, one may be interested in evaluating 
various functions f(Ai, . . . , Ar) of them, which may involve fit procedures and other 
complicated operations. The standard stochastic estimate of such an observable is 

f = f(a 1 ,...,a R ) (A.5) 
and the associated series of jackknife estimates is defined by 

// = /«i,...,4,i), i = l,...,N. (A.6) 
A little algebra then shows that the expression 

^ 2 = E(//-/) 2 (A.7) 

i=i 

provides an estimate of the statistical variance of /, which coincides with the usual 
error propagation formula (the one that involves the gradient of /) up to terms of 
order 1/N. Similarly the error covariance of / and any other function g is obtained 
by summing (// — f){gf — g) over the jackknife samples. 

In practice the error formula (A.7) proves to be very convenient. If an observable is 
a function of previously calculated observables, for example, one can take advantage 
of the fact that the composition of functions is associative, i.e. the jackknife series // 
is simply obtained by inserting the jackknife series of the arguments, independently 
of whether these are primary or not. The data analysis can thus proceed in steps, 
starting from the primary observables and progressing to more and more complicated 
observables. 
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A . 3 Combining data from different runs 

Simulations of lattice QCD at different sea-quark masses, lattice spacings, etc., can 
be assumed to be statistically independent. The statistical variance of any observable 
that depends on data from several simulations is therefore the sum of the associated 
partial variances. This rule can easily be accommodated in the jackknife analysis 
by embedding the jackknife series of the observables in extended series that include 
all simulations on which the observable depends. 

The method is best explained by considering two simulations, where N\ measure- 
ments of some observables A r are made in the first and N2 measurements of some 
other observables B s in the second. The associated jackknife series a^ l5 . . . ,a^ Nl 
and b J s 1 , . . . , b J s N are then computed as before, starting from the primary observ- 
ables in each simulation. Next they are embedded in extended series 

a J r ^,...,a^ Ni , a r ,...,a r ^ and b s , . .. , b s b J sl , . . . , b J sN2 (A.8) 

N2 elements Ni elements 

of length Ni + N2 such that the first N\ elements are occupied by the jackknife 
series from the first simulation and the last N2 elements by those from the second 
simulation. 

With this assignment, and if the extended series are treated as ordinary jackknife 
series, the correct error correlation matrix of the full set Ai, . . . , Ar, Bi, ... , Bs of 
observables is obtained. Moreover, we may define the jackknife series of any ob- 
servable f(Ai, . . . Ar, Bi, ... , Bs) in the standard manner and compute its variance 
using eq. (A. 7). The embedding trick thus allows the statistical errors to be propa- 
gated as if there were a single simulation. 



Appendix B. Properties of the auxiliary function h{t;E,E') 

The symmetry properties 

h(t; E, E') = h{T - t; E, E') = h(t; E' , E) (B.l) 

are an immediate consequence of the definition (5.9) of the function h(t;E,E'). It 
is also straightforward to verify that 

h(t;E,E') = 2e-^ B+i5 ') T cosh (\{E' - E)(T - 2t)) , (B.2) 
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and h(t; E, E') is thus a convex function of t which attains its minimum at t = \T. 

B.l Product inequality 

We now show that the inequality 

h(t; E, E' + E") < h(t; E, E')h(t; 0, E") (B.3) 
holds for all values of the arguments t, E, E' and E" . To this end, first note that 

cosh(a + 0) < cosh(a + 0) + cosh(a — 0) =2 cosh a cosh 0. (B.4) 
Substituting a = \ {E' - E)(T - 2t) and = \E"(T - 2t), this inequality becomes 

cosh (\{E' + E" -E)(T- 2t)) < 

2 cosh (\{E' - E)(T - 2t)) cosh {\E"{T - 2t)) , (B.5) 
which is easily seen to coincide with (B.3) after inserting the representation (B.2). 
B.2 Monotonicity property 

If t and s are in the range s < t < |T, and if M > 0, it follows from eq. (B.2) that 
the ratio 

h(s;0,M) 



h(t;0,M) 



(B.6) 



is greater than 1. A less obvious statement is that the ratio increases monotonically 
from r = 1 to r = oo when M goes from zero to infinity. 

In order to show this, we insert eq. (B.2) and work out the quotient 

q = Lzi = tanh (\M{t - a)) tanh {\M(T - t - s)) . (B.7) 

In the specified range of t and s, the arguments of the hyperbolic functions in this 
equation are non-negative and monotonically growing with M. The quotient thus 
rises monotonically from to 1 when M goes from zero to infinity, which proves our 
claim, since r and q are monotonically related to each other. 
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Appendix C. Tables of meson masses and decay constants 

The simulation results tabulated in this appendix were obtained following the lines 
of sect. 5. In all cases, the r quark was taken to be a sea quark, i.e. the associated 
hopping parameter n r was set to K sea . The hopping parameter n s of the other quark, 
on the other hand, ranged over 4 or 5 values, one of which being n sea _. 

For each series of runs, we quote the quark mass sums m rs , the pseudo-scalar 
meson masses Mps and matrix elements Gps, and the vector meson masses My, all 
given in lattice units (tables 5, 7 and 9). Some combinations of these quantities are 
printed in tables 6, 8 and 10. The errors given in brackets are statistical only. 

If so desired, the quoted results can be converted to physical units by substituting 
the estimates 0.0717(15), 0.0521(7) and 0.0784(10) fm for the spacings of the A, B 
and D lattices [9]. The quark mass sums m rs , the matrix elements Gps and the 
decay constants Fps then also need to be renormalized, 

m rs — ► Z^Z^rrirs, Gps — > ZpGps, F PS — > Za-Fps, (CI) 

where Za and Zp denote the (mass-independent) renormalization constants of the 
non-singlet axial current and density. Moreover, in order to guarantee the O(a) im- 
provement of these quantities in the improved theory, the renormalization constants 
must be modified according to 

Z x -> Z x (l + b x am sea + \b x am rs ) , (C.2) 

with properly adjusted coefficients 6x and b x [17,33] (the figures quoted in tables 9 
and 10 include the contribution of the operator improvement term proportional to 
ca but not the 1 + 0(am) renormalization factors). 
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Table 5. Results for m rs , Mps, Gps and My (lattices A\ — A3) 



Run 




K s 


am rs 


aM PS 


a 2 G PS 


aM v 


M 


0.15750 


0.15750 
0.15800 
0.15825 
0.15835 


0.0548(5) 
0.0472(6) 
0.0434(6) 
0.0419(6) 


0.2726(19) 
0.2536(19) 
0.2438(20) 
0.2398(21) 


0.0881(12) 
0.0859(12) 
0.0848(13) 
0.0844(13) 


0.389(4) 
0.379(5) 
0.373(5) 
0.371(5) 


A 2 


0.15800 


0.15750 
0.15800 
0.15825 
0.15835 


0.0359(3) 
0.0285(3) 
0.0249(3) 
0.0235(3) 


0.2137(18) 
0.1913(19) 
0.1790(21) 
0.1738(22) 


0.0703(12) 
0.0682(13) 
0.0671(14) 
0.0666(15) 


0.344(3) 
0.334(4) 
0.329(5) 
0.328(5) 


A 3 


0.15825 


0.15750 
0.15800 
0.15825 
0.15835 


0.0281(4) 
0.0208(4) 
0.0172(4) 
0.0158(4) 


0.185(3) 
0.160(3) 
0.147(4) 
0.141(4) 


0.0617(19) 
0.0599(22) 
0.0593(23) 
0.0592(24) 


0.327(5) 
0.317(7) 
0.312(8) 
0.311(9) 
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Table 6. Combinations of m rs , Mps, Gps and My (lattices A\ — A3) 



Run 






aMp S /m rs 


aFp S 


Fps/Mv 


M 


0.15750 


0.15750 


1.357(17) 


0.0650(7) 


0.1669(21) 






0.15800 


1.363(19) 


0.0630(7) 


0.1664(24) 






0.15825 


1.369(21) 


0.0619(8) 


0.166(3) 






0.15835 


1.372(22) 


0.0615(8) 


0.166(3) 


A 2 


0.15800 


0.15750 


1.272(20) 


0.0553(7) 


0.161(3) 






0.15800 


1.282(25) 


0.0532(7) 


0.159(3) 






0.15825 


1.29(3) 


0.0522(8) 


0.159(3) 






0.15835 


1.29(3) 


0.0518(8) 


0.158(4) 




0.15825 


0.15750 


1.22(4) 


0.0505(8) 


0.154(3) 






0.15800 


1.23(5) 


0.0486(10) 


0.153(4) 






0.15825 


1.25(6) 


0.0474(11) 


0.152(5) 






0.15835 


1.26(7) 


0.0469(12) 


0.151(6) 
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Table 7. Results for m rs , Mps, Gps and My (lattices B\ — B 4 ) 



Run 




K s 


am rs 


aM PS 


a 2 G PS 


aM v 


B 1 


0.15410 


0.15410 
0.15425 
0.15440 
0.15455 


0.03889(18) 
0.03631(18) 
0.03375(18) 
0.03120(18) 


0.1958(9) 
0.1892(9) 
0.1824(10) 
0.1754(10) 


0.0453(5) 
0.0447(5) 
0.0441(6) 
0.0435(6) 


0.2896(17) 
0.2858(17) 
0.2820(18) 
0.2782(19) 


B 2 


0.15440 


0.15410 
0.15425 
0.15440 
0.15455 


0.02696(13) 
0.02440(14) 
0.02187(14) 
0.01935(14) 


0.1619(11) 
0.1546(12) 
0.1470(12) 
0.1391(13) 


0.0384(7) 
0.0379(7) 
0.0374(7) 
0.0370(8) 


0.2518(21) 
0.2475(22) 
0.2432(24) 
0.239(3) 


B 3 


0.15455 


0.15410 
0.15425 
0.15440 
0.15455 


0.02185(12) 
0.01927(12) 
0.01668(12) 
0.01409(13) 


0.1416(12) 
0.1329(13) 
0.1235(14) 
0.1132(15) 


0.0333(7) 
0.0326(7) 
0.0318(7) 
0.0310(8) 


0.2418(24) 
0.238(3) 
0.233(3) 
0.230(3) 




0.15462 


0.15410 
0.15425 
0.15440 
0.15455 
0.15462 


0.02029(16) 
0.01774(16) 
0.01521(17) 
0.01269(17) 
0.01151(17) 


0.1328(10) 
0.1242(11) 
0.1151(12) 
0.1055(14) 
0.1008(15) 


0.0317(6) 
0.0312(6) 
0.0307(6) 
0.0302(7) 
0.0300(8) 


0.237(3) 
0.233(3) 
0.229(3) 
0.224(4) 
0.223(4) 
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Table 8. Combinations of m rs , Mps, Gps and My (lattices B\ — £> 4 ) 



Run 






aMp S /m rs 


aFps 


Fps/Mv 


B 1 


0.15410 


0.15410 
0.15425 
0.15440 
0.15455 


0.986(7) 
0.986(8) 
0.986(8) 
0.986(9) 


0.0460(4) 
0.0453(4) 
0.0447(4) 
0.0441(5) 


0.1587(17) 
0.1586(18) 
0.1585(19) 
0.1584(20) 


B 2 


0.15440 


0.15410 
0.15425 
0.15440 
0.15455 


0.973(14) 
0.979(15) 
0.988(17) 
1.000(19) 


0.0395(4) 
0.0387(4) 
0.0379(4) 
0.0370(4) 


0.1567(19) 
0.1562(20) 
0.1556(21) 
0.1549(23) 


B 3 


0.15455 


0.15410 
0.15425 
0.15440 
0.15455 


0.918(15) 
0.916(17) 
0.914(19) 
0.910(22) 


0.0363(3) 
0.0356(4) 
0.0348(4) 
0.0340(4) 


0.1502(21) 
0.1497(23) 
0.149(3) 
0.148(3) 




0.15462 


0.15410 
0.15425 
0.15440 
0.15455 
0.15462 


0.869(13) 
0.869(15) 
0.871(18) 
0.878(23) 
0.88(3) 


0.0365(4) 
0.0359(4) 
0.0352(5) 
0.0344(5) 
0.0340(6) 


0.154(3) 
0.154(3) 
0.154(3) 
0.153(4) 
0.153(4) 
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Table 9. Results for m rs , Mps, Gps and My (lattices D\ — D 5 ) 



Run 




K s 


am rs 


aM PS 


a 2 G PS 






0.13550 


0.13550 


0.06771(21) 


0.3286(10) 


0.1069(15) 


0.464(3) 






0.13590 


0.05704(22) 


0.3017(10) 


0.1030(15) 


0.447(3) 






0.13610 


0.05165(23) 


0.2873(11) 


0.1008(15) 


0.438(3) 






0.13620 


0.04893(24) 


0.2799(12) 


0.0998(15) 


0.434(4) 


D 2 


0.13590 


0.13550 


0.04968(13) 


0.2758(8) 


0.0920(11) 


0.4173(24) 






0.13590 


0.03914(14) 


0.2461(9) 


0.0891(11) 


0.401(3) 






0.13610 


0.03383(14) 


0.2301(9) 


0.0880(12) 


0.394(4) 






0.13620 


0.03112(15) 


0.2218(10) 


0.0878(13) 


0.390(4) 


D 3 


0.13610 


0.13550 


0.04092(14) 


0.2440(10) 


0.0811(12) 


0.382(3) 






0.13590 


0.03041(14) 


0.2110(11) 


0.0780(13) 


0.363(4) 






0.13610 


0.02514(15) 


0.1929(12) 


0.0766(14) 


0.354(5) 






0.13620 


0.02249(15) 


0.1832(13) 


0.0760(15) 


0.349(5) 


D 4 


0.13620 


0.13550 


0.03728(14) 


0.2335(11) 


0.0813(11) 


0.374(4) 






0.13590 


0.02686(15) 


0.1993(12) 


0.0785(12) 


0.356(4) 






0.13610 


0.02168(15) 


0.1800(13) 


0.0771(13) 


0.348(5) 






0.13620 


0.01909(15) 


0.1695(14) 


0.0765(13) 


0.345(6) 




0.13625 


0.13550 


0.03474(13) 


0.2249(11) 


0.0784(13) 


0.376(5) 






0.13590 


0.02428(13) 


0.1881(11) 


0.0747(14) 


0.359(6) 






0.13610 


0.01910(13) 


0.1672(13) 


0.0729(15) 


0.350(7) 






0.13620 


0.01651(14) 


0.1559(14) 


0.0722(16) 


0.346(8) 






0.13625 


0.01522(14) 


0.1499(15) 


0.0719(17) 


0.344(9) 
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Table 10. Combinations of m rs , Mps, Gps and My (lattices D\ — D 5 ) 



Run 






aMp S /m rs 


aFp S 


Fps/Mv 


D 1 


0.13550 


0.13550 
0.13590 
0.13610 
0.13620 


1.594(9) 
1.596(10) 
1.598(11) 
1.601(12) 


0.0671(9) 
0.0645(9) 
0.0631(9) 
0.0624(9) 


0.1445(20) 
0.1444(21) 
0.1440(23) 
0.1438(23) 


D 2 


0.13590 


0.13550 
0.13590 
0.13610 
0.13620 


1.531(9) 
1.547(10) 
1.565(12) 
1.581(14) 


0.0601(6) 
0.0576(7) 
0.0562(7) 
0.0556(7) 


0.1441(17) 
0.1435(20) 
0.1428(22) 
0.1424(24) 


D 3 


0.13610 


0.13550 
0.13590 
0.13610 
0.13620 


1.454(11) 
1.465(14) 
1.480(17) 
1.492(19) 


0.0558(6) 
0.0533(7) 
0.0518(7) 
0.0510(8) 


0.1461(20) 
0.1467(23) 
0.146(3) 
0.146(3) 


D 4 


0.13620 


0.13550 
0.13590 
0.13610 
0.13620 


1.462(13) 
1.478(17) 
1.494(20) 
1.505(23) 


0.0556(6) 
0.0531(7) 
0.0516(7) 
0.0508(7) 


0.1487(22) 
0.149(3) 
0.148(3) 
0.147(3) 


D 5 


0.13625 


0.13550 
0.13590 
0.13610 
0.13620 
0.13625 


1.456(15) 

1.457(19) 

1.464(24) 

1.47(3) 

1.48(3) 


0.0539(8) 
0.0512(8) 
0.0498(8) 
0.0490(9) 
0.0487(9) 


0.143(3) 
0.143(3) 
0.142(3) 
0.142(4) 
0.142(4) 
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